# load data
load("[INSERT PATH]/EMN_replication_JEPP.RData")

### load library
library(statnet)
library(ggplot2)
library(tidyverse)
library(stargazer)
library(gridExtra)

#######################################################################################
################################## Figures and Tables ################################## 

### Figure 1 ###
p1 <- ggplot(data=queries_trend, aes(x=as.factor(Year))) +
  geom_bar(aes(y=Queries), stat="identity", position="identity", alpha=.8, fill='black', color='black') +
  xlab("Year") +
  ylab("Number of queries") +
  theme_minimal()

p2 <- ggplot(data=queries_trend, aes(x=as.factor(Year))) +
  geom_bar(aes(y=Interactions), stat="identity", position ="identity", alpha=.3, fill='grey10', color='grey10') +
  xlab("Year") +
  ylab("Number of answers") +
  theme_minimal()

grid.arrange(p1, p2, nrow = 1)

### Figure 1 appendix ###
ggplot(data= queries_theme) + 
  geom_bar(aes(x=Theme, y=Queries),stat="identity") +
  coord_flip() +
  theme_minimal()
dev.off()

### Figure 2 appendix ###
p3 <- ggplot(data= queries_by_ms) + 
  geom_bar(aes(x = reorder(MS, Queries), y=Queries),stat="identity") +
  xlab("Member State") +
  ylab("Number of questions") +
  coord_flip() +
  theme_minimal()

p4 <- ggplot(data= queries_by_ms) + 
  geom_bar(aes(x = reorder(MS, Answers), y=Answers),stat="identity") +
  xlab("Member State") +
  ylab("Number of answers") +
  coord_flip() +
  theme_minimal()

grid.arrange(p3, p4, nrow = 1)
dev.off()


### Table 2 in appendix ###
stargazer(df,  type="text", median=TRUE)

################################## ERGM ##################################### 
model <- ergm(full_network ~ sum +
                transitiveweights("min", "max", "min") + 
                cyclicalweights("min", "max", "min") +
                mutual("min") +
                #nodecov("gov_restrictiveness") + #H1
                absdiff("admin_capacity_2015_2018") + #H2
                #nodefactor("affected") + #H3
                nodecov("asylum_applications") +
                nodefactor("affected2") +
                nodecov("tolerance_neg_change") + #H4
                nodecov("NCP_size") +
                nodematch("NCP_model", diff=FALSE) +
                edgecov("2016_dummy", attrname="2016 queries (dummy") +
                edgecov("2017_dummy", attrname="2017 queries (dummy"),
              response = "relations", reference = ~Binomial(11))

### Table 1 ###
stargazer(model,  type="text", single.row=FALSE)

################################## ERGM ##################################### 
model2 <- ergm(full_network ~ sum +
                transitiveweights("min", "max", "min") + 
                cyclicalweights("min", "max", "min") +
                mutual("min") +
                #nodecov("gov_restrictiveness") + #H1
                #absdiff("admin_capacity_2015_2018") + #H2
                nodematch("admin_capacity_dummy", diff=TRUE) +
                #nodefactor("affected") + #H3
                nodecov("asylum_applications") +
                nodefactor("affected2") +
                nodecov("tolerance_neg_change") + #H4
                nodecov("NCP_size") +
                nodematch("NCP_model", diff=FALSE) +
                edgecov("2016_dummy", attrname="2016 queries (dummy") +
                edgecov("2017_dummy", attrname="2017 queries (dummy"),
              response = "relations", reference = ~Binomial(11))

### Table 2 ###
stargazer(model,model2,  type="text", single.row=FALSE)

########################## mcmc diagnostics #############################

par(mar=c(3,3,3,3))

mcmc.diagnostics(model,vars.per.page = 6)
dev.off()


######################### multicolinearity check ########################

VIF.ERGM(model) #see function below
VIF.ERGM(model2) #see function below


######################### interpret coefficients #########################

exp(coef(model)[5]) #h2: effectiveness similarity
exp(coef(model)[6]) #h3: problem pressure - nr of applications
exp(3.728241*coef(model)[6]) - exp(602.146354*coef(model)[6]) #change from lowest to highest range
exp(coef(model)[7]) #h3: problem pressure - affectedness - destination states
exp(coef(model)[8]) #h3: problem pressure - affectedness - first-entry state (not significant)
exp(coef(model)[9]) #h4: change in tolerance
exp(-7.586 *coef(model)[9]) - exp(13.649*coef(model)[9]) #change from lowest to highest range
exp(coef(model)[10]) # capacity
exp(coef(model)[12]) # 2016
exp(coef(model)[13]) # 2017

########################### goodness of fit ############################

## Simulate from model fit:
modelv.sim <- simulate(model,
                       nsim = 10000, stats=TRUE)

## Create density plots for each parameter
par(mfrow = c(7,2),mar=c(2, 2, 2, 2) +0.1)

 # density
vnet.obs <- model$target.stats[1]
plot(density(modelv.sim[,1]), main="Sum of edges")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 1])
min(mean(modelv.sim[, 1]>vnet.obs) ,mean(modelv.sim[, 1]<vnet.obs))*2

# transitivity
vnet.obs <- model$target.stats[2]
plot(density(modelv.sim[,2]), main="Transitivity")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 2])
min(mean(modelv.sim[, 2]>vnet.obs) ,mean(modelv.sim[, 2]<vnet.obs))*2

# cyclical weights
vnet.obs <- model$target.stats[3]
plot(density(modelv.sim[,3]), main="Cyclical triads")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 3])
min(mean(modelv.sim[,3]>vnet.obs) ,mean(modelv.sim[, 3]<vnet.obs))*2

# reciprocity
vnet.obs <- model$target.stats[4]
plot(density(modelv.sim[,4]), main="Reciprocity")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 4])
min(mean(modelv.sim[, 4]>vnet.obs) ,mean(modelv.sim[,4]<vnet.obs))*2


# Government effectiveness
vnet.obs <- model$target.stats[5]
plot(density(modelv.sim[,5]), main="Government effectiveness dissimilarity ")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 5])
min(mean(modelv.sim[, 5]>vnet.obs) ,mean(modelv.sim[, 5]<vnet.obs))*2

# problem pressure - nr of applicants
vnet.obs <- model$target.stats[6]
plot(density(modelv.sim[,6]), main="Number of applicants")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 6])
min(mean(modelv.sim[, 6]>vnet.obs) ,mean(modelv.sim[, 6]<vnet.obs))*2

# problem pressure - destination states
vnet.obs <- model$target.stats[7]
plot(density(modelv.sim[,7]), main="Destination states")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 7])
min(mean(modelv.sim[, 7]>vnet.obs) ,mean(modelv.sim[, 7]<vnet.obs))*2

# problem pressure - first-entry states
vnet.obs <- model$target.stats[8]
plot(density(modelv.sim[,8]), main="First-entry states")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 8])
min(mean(modelv.sim[, 8]>vnet.obs) ,mean(modelv.sim[, 8]<vnet.obs))*2

# Decrease in immigration tolerance
vnet.obs <- model$target.stats[9]
plot(density(modelv.sim[,9]), main="Decrease in immigration tolerance")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 9])
min(mean(modelv.sim[, 9]>vnet.obs) ,mean(modelv.sim[9]<vnet.obs))*2

# NCP personnel capacity
vnet.obs <- model$target.stats[10]
plot(density(modelv.sim[,10]), main="NCP personnel capacity")
abline(v=vnet.obs)
density.default(x = modelv.sim[,10])
min(mean(modelv.sim[, 10]>vnet.obs) ,mean(modelv.sim[, 10]<vnet.obs))*2

# NCP
vnet.obs <- model$target.stats[11]
plot(density(modelv.sim[,11]), main="Similar NCP")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 11])
min(mean(modelv.sim[, 11]>vnet.obs) ,mean(modelv.sim[, 11]<vnet.obs))*2

# Dummy 2016
vnet.obs <- model$target.stats[12]
plot(density(modelv.sim[,12]), main="Year 2016")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 12])
min(mean(modelv.sim[, 12]>vnet.obs) ,mean(modelv.sim[, 12]<vnet.obs))*2

# Dummy 2017
vnet.obs <- model$target.stats[13]
plot(density(modelv.sim[,13]), main="Year 2017")
abline(v=vnet.obs)
density.default(x = modelv.sim[, 13])
min(mean(modelv.sim[, 13]>vnet.obs) ,mean(modelv.sim[, 13]<vnet.obs))*2

dev.off()

######################### VIF function #########################

VIF.ERGM<-function(my.ergm){
  require(ergm)
  
  cor.mat<-cov2cor(my.ergm$covar) #calculate correlation matrix
  corr5<-cor.mat[-c(1),-c(1)] ##omit edges term
  
  corr5<-corr5[!is.na(corr5[1:nrow(corr5)]),]
  corr5<-corr5[,which(!is.na(corr5[1,1:ncol(corr5)]))]
  
  VIFS<-matrix(0,nr=1,nc=ncol(corr5))
  
  for(i in 1:ncol(corr5)){
    
    gvec<-as.vector(corr5[-c(i),i]) ##create vector of correlations between covariate of interest and other covariates in the model
    tgvec<-t(gvec)            
    xcor<-solve(corr5[-c(i),-c(i)]) ##create square matrix of correlations between covariates in the model other than the one of interest
    Rsq<-tgvec%*%xcor%*%gvec
    VIFS[1,i]<-1/(1-Rsq)
  }
  ##Columns are covariates as they appear in the ERGM object
  colnames(VIFS)<-names(my.ergm$coef[-c(1)])
  message("Higher values indicate greater correlation.\nVIF > 20 is concerning, VIF > 100 indicates severe collinearity.")
  VIFS
}
